Lattice Boltzmann Method with regularized non-equilibrium distribution functions 
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A new lattice Boltzmann (LB) model is introduced, based on a regularization of the pre-collision 
distribution functions in terms of the local density, velocity, and momentum flux tensor. The model 
dramatically improves the precision and numerical stability for the simulation of fluid flows by LB 
methods. This claim is supported by simulation results of some 2D and 3D flows. 

PACS numbers: 47.11. +j, 05.20.Dd 



The lattice Boltzmann (LB) model is a recent tech- 
nique for the simulation and modeling of fluid flows [1- 
5], During the past fifteen years it has been successfully 
applied to many challenging problems in hydrodynamics 
as well as reaction-diffusion processes and wave propa- 
gation phenomena [6] . A particular shortcoming of this 
technique are numerical instabilities, that may develop 
at high Reynolds numbers. Several improvements to the 
method have been proposed, which however either induce 
a substantial complication of the original algorithm, or 
require a cumbersome fine-tuning of adjustable parame- 
ters [7, 8]. The present paper introduces a new method 
which fits quite naturally into the framework of classi- 
cal LB models and offers both increased accuracy and 
stability at very low cost. 

The LB approach considers a mesoscopic description 
of the fluid on a regular lattice of spacing Sr in d- 
dimensions. The central quantities of the LB approach 
are distribution functions fi(r,t), which denote the den- 
sity of particles entering a lattice site r at discrete time t 
with velocity Vi. The Vi are vectors connecting any lattice 
site r with its z neighbors f+Stvi, St being the time step 
and z the lattice coordination number. A vector vq = 
corresponding to a rest population /o is also introduced. 
The LB dynamics are expressed as 

Mr + Stv h t + St) = fi(f, t) + fi<(/(j* t)), (1) 

where i, here and in subsequent formulas, runs from to 
z. The dynamics can be split conceptually into a collision 
step by defining f° ut = fi + fij(/), and a propagation 
step: fi(f+ 5tvi,t + St) — f° ut (r,t). During the collision 
step, the advected particle streams fi are summed up 
with the collision terms f2i, which are given functions of 
the fi's. They describe how fluid particles colliding at site 
r change their velocities to vi. Then, at the propagation 
step, the fluid particles are streamed to the neighboring 
site f+StHi. 

As in any standard kinetic theory, the macroscopic 
quantities are obtained by taking the first velocity mo- 
ments of the distribution functions: 



where p, u, and II are the fluid density, momentum, and 
momentum flux tensor respectively (Note that the ac- 
tual momentum flux tensor in LB models has an extra 
lattice contribution, which adds on to II.). Here and in 
what follows, Greek indices label the components of two- 
dimensional (2D) resp. three-dimensional (3D) physical 
space, whereas Latin indices refer to the z+l-dimensional 
space of the distribution functions. Vectors situated in 
the former space are characterized by an arrow on top 
of the letter, and in the latter space, simply by omitting 
the index. 

The collision term O is chosen in such a manner that 
mass and momentum are conserved exactly (without dis- 
cretization error), so as to closely reflect the physical laws 
at the base of hydrodynamics. Its most common imple- 
mentation, the BGK model, expresses a single-time relax- 
ation to a given local equilibrium function f eq , depend- 
ing only on the conserved quantities p and u calculated 
from (2): 



(3) 



where < u) < 2 is the relaxation parameter, directly 
related to the dynamic fluid viscosity v. 

The expression for f eq comes from a low Mach num- 
ber truncated Maxwell-Boltzmann distribution and is ad- 
justed to obtain the correct momentum flux tensor: 



n «/3 = X! fk^VkaVip = PC 2 s S a p + pU a U/3, (4) 

where c s is the speed of sound. The equilibrium term f eq 
reads [1-3] 
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(5) 



PU = ^2 fk^ki n Q/3 = 22 fkVkaVkp, (2) 
k=0 k=0 



where a repeated Greek index implies a summation over 
this index. The tensors Qi a p are defined to be Qi a p = 
ViaVi/3 — c^Sap, and the ij's, as well as c s , are coefficients 
specific to the lattice topology. 

The connection between the LB method and the cor- 
responding hydrodynamics is obtained through a Tay- 
lor expansion, up to second order in St, of the fi- 
nite differences in the left hand side of Eq. (1), and a 
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multiscale Chapman-Enskog expansion / = YlkLaf^' 
The zeroth-order term yields the equilibrium distribution 
value /(°) = f eq , and the remaining terms are denoted 
as f neq : 

fi eq = f i -f? t and n ne9 = n-n e<? . (6) 

For the BGK model, the first-order multiscale Chapman- 
Enskog procedure gives [3] 

~ /i^ = %tiQia/3d a pU0, (7) 

UJC S 

and we obtain 

z 2 

U a7 ~ ^2 f^^kaVkP = " {daPUp + dppU a ) . (8) 

H — LO 

k=0 

Using expressions (4) and (8) together with the lattice 
contribution to the momentum flux (see for instance [3]), 
it can be shown that u obeys the Navier-Stokes equation 
with the viscosity given by 

However, in actual numerical simulations, the pro- 
posed theoretical description of the LB dynamics is not 
fully obeyed because St and Sr are not arbitrarily small, 
and also because higher order derivatives are neglected 
in the approximation (7). As a result, the numerical 
behavior departs from its hydrodynamic limits and nu- 
merical instabilities may appear if some quantities vary 
too sharply over time and space. 

The inaccuracy of the first-order terms /W becomes 
apparent, e.g., upon the observation that, according to 
Eq. (7), is symmetric with respect to spatial reflec- 
tions: the difference — /j 1 ' vanishes along directions 
i,j for which = —Vj. In practice, this relation is not 
necessarily obeyed by the non-equilibrium parts of the 
distribution functions. On Fig. 1, Ay = /" e9 — /" e<? is 
plotted, for a given couple {i,j}, on ground of some nu- 
merical data of the Kovasznay flow described below. It 
appears to take nonnegligible values at the scale of non- 
equilibrium terms (up to 30%). 

To reduce the discrepancy between f neq and f^ 1 ', we 
propose a regularization procedure whose goal is to force 
the numerical scheme to comply as much as possible with 
the theoretical framework exposed above. For this pur- 
pose we recompute f neq prior to the collision step so as 
to enforce f neq = The key of our regularization 

procedure is the observation that Eqs. (7) and (8) can be 
combined to give 

In conclusion, our regularization procedure amounts to 
computing the regularized values j' 1 ' of f neq according 




FIG. 1: Difference between non-equilibrium parts of distri- 
bution functions along two opposite directions vi/v = (1,0) 
and Vz/v = ( — 1,0). Each data point is obtained from one 
lattice site on a numerical simulation of a Kovasznay flow, at 
Re = I. 

to the following steps: 

I Eq ^ ] (^)^ 5) / e «M (ID 

Eq. (6) / f neq \ Eq. (10) (1) 

' ^n ne « ) * 1 

Then, the standard BGK collision is applied to / = f eq + 
/W, and the regularized collision step of the dynamics 
reads 

fr t =f?+o--x>)ti 1) . (12) 

Note that since J^k Qkap = J2k Qkapvk = 0, the above 
scheme still conserves mass and momentum exactly . 

In order to better understand the way the steps de- 
scribed in Eq. (11) act on the distribution functions, it 
is illuminating to study the dynamics in the (z + 1)- 
dimensional space of the velocity moments. These mo- 
ments [a few of them are shown in Eq. (2)] are associ- 
ated in kinetic theory with so-called modes of the colli- 
sion operator and can be related to transport phenom- 
ena during the collision process. In general, the moment 
space is related to the space of the distribution functions 
through an invertible linear mapping whose matrix M 
is explicited, e.g., in [7]. The regularized dynamics pre- 
sented in Eqs. (11, 12) can be reformulated as 

f out = f eq + {l-Lo)M- l AMf neq , (13) 

where A = MRM , with i?y = ^tsQiapCjaCjp. In 
2D and under the assumption of fluid incompressibil- 
ity Tr(II) = 0, the matrix A is found to be diagonal: 
A = diag(\o, Xi, ■ ■ ■ , X z ), where = A,- 2 = 1 for the 
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components of the momentum flux tensor[17] and A^ = 
for the other moments. In the general (compressible) 
case, additional off-diagonal contributions appear in the 
energy and square-energy moments. This interpretation 
of the dynamics shows that, except for compressibility 
effects, the regularized dynamics directly kills all modes 
but the ones associated to the momentum flux tensor. 

It is interesting to compare the regularized method 
with so-called multi-relaxation-time (MRT) models, [7, 
9], which propose the following general formulation of 
the LB dynamics: f out = f(r,t) - M^S^T - T eq ), 
where T = Mf is the moment space representation 
of the distribution functions, S is a diagonal matrix 
S = diag{sQ, S\, ■ • ■ , s z ) containing z + l individual relax- 
ation parameters s,-, and T eq , the equilibrium distribu- 
tion in moment space, depends on a set of adjustable pa- 
rameters. By fixing those adjustable parameters through 
the relation T eq — Mf eq and the relaxation parameters 
through Si — lo for all non-conservative momenta, the 
usual BGK dynamics are recovered. It has however been 
argued [7], that the stability of the BGK scheme is en- 
hanced by an appropriate choice of the various relaxation 
parameters. 

In the case T eq = Mf eq , the MRT model takes the 
following form: f out = f(>r,t)-M- 1 SM(f-f e ''), which, 
in analogy with Eq. (13), can be reformulated as 

j-out = jeq + (J. _ M^SM)/^. (14) 

Here, the identity term 1 is due to the advected dis- 
tribution functions, which are not touched upon by the 
MRT correction to the BGK model. Eqs. (13, 14) make 
the main difference between our regularized model and 
the MRT approach apparent: while in the MRT ap- 
proach [Eq. (14)] non-physical modes are relaxed to a 
local equilibrium inside the collision term, in the regu- 
larized model [Eq. (13)] these modes are more radically 
eliminated in both the advected particles and the colli- 
sion term. Therefore, when increasing the stability in the 
simulation of a Navier-Stokes fluid How is the only issue, 
our method is comparatively simpler from a theoretical 
viewpoint, and efficient to implement. 

We now turn to numerical verifications of the regular- 
ized model on two 2D flows using a D2Q9 lattice, and one 
3D flow using a D3Q19 lattice [1-3]. The first test con- 
cerns the simulation of a Kovasznay flow, which approx- 
imates the stationary 2D flow behind a regular grid. An 
analytical solution for this flow, proposed in [10], takes 
the following form: 

u x = itoo(l — exp(Ax/-L) • cos(27rj//L)) and (15) 

Uy — Woo — exp(Ax/L) ■ sin(27ry/i), with 

A = Re/2 - v/47r 2 + i?e 2 /4, 

where Uoo is the asymptotic velocity of the fluid, Re = 
UqoL/v is the Reynolds number, and L defines the length 
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FIG. 2: Relative error of numerical result on Kovasznay flow 
in wake region. Both traditional BGK and the regularized 
model are tested on two common boundary conditions. 

scale of the problem. The simulations are performed in 
the wake of the grid, in the intervals x G [L/2,2L] and 
y G [-L/2, 3L/2], with Re = 10, = 0.01 v, and with a 
varying grid resolution N — L/Sr. Keeping the velocity 
constant in terms of the lattice unit v amounts to fixing 
the Mach number Ma — u/c s at a value sufficiently small 
to mimic an incompressible flow. Given that the flow is 
periodic in y-direction, the upper and the lower bound- 
ary of the simulation can be chosen periodic, whereas the 
Kovasznay solution [Eq. (15)] is imposed through Dirich- 
let boundary conditions on the left and right boundary. 
After the simulation has stabilized, the numerical result 
is compared with the solution [Eq. (15)] through an Li 
norm on each grid point, and then averaged over space. 
The result is shown in Fig. 2, on two commonly used 
implementations of the boundary conditions (be); be (1) 
[11] and be (2) [12]. The accuracy of the simulation with 
respect to the grid resolution is of order 2 to 2.5 when the 
BGK model is used, whereas the regularized model is al- 
most third-order accurate. On the BGK simulations with 
be (1), data points for small grids are missing because nu- 
merical instabilities make them impossible, whereas the 
regularized model has no such stability deficiencies. 

The second test case implements a flow in a 2D square 
cavity whose top-wall moves with a uniform velocity. 
Both standard BGK and the regularized model are first 
compared with the reference solution of Ghia e.a. [13], 
on a lattice size of N x N with N = 129, at Re = 100 
and a top- wall velocity uq = 0.02 v. A boundary condi- 
tion described in [14] is used. The reference solution [13] 
proposes a set of accurate numerical values for some x- 
and some y-components of the velocity on chosen space 
points. An L\ norm error with respect to these reference 
points is averaged over all available points and normal- 
ized with respect to uq. For the BGK model, this yields 



4 



1000 




50 60 70 80 90 100 110 120 130 
Grid resolution N 



FIG. 3: Simulation of 2D cavity flow for fixed Mach number, 
o, *: maximal stable Reynolds number, numerically deter- 
mined; solid line: least-square linear fit of the data points 
(parameters of the fit are indicated on the graph). 



numerically verified to fit the predictions of the theory 
of fluid turbulence. However, at a smaller Kolmogorov 
length (and thus higher Re) Ik = 0.06 5r, BGK is numeri- 
cally unstable, whereas numerical stability is still ensured 
by the renormalized model. This observation suggests 
that the physics of the small scales are represented more 
accurately by the renormalized model than by BGK. 

In this paper, a novel numerical scheme has been pre- 
sented for the simulation of fluid flows by the LB method. 
It has been compared with the traditional BGK method 
and shown to be substantially more precise on a prob- 
lem with mathematically well defined boundaries, dra- 
matically more stable on a problem with high pressure 
gradients on a critical point, and more robust against 
an excessive energy input in a turbulent flow. Given its 
conceptual simplicity, we highly recommend its use as 
an alternative model for the simulation of complex fluid 
flows. We thankfully acknowledge the support by the 
Swiss National Science Foundation (SNF). 



an error of e = 3.71 TO -3 , and for the regularized method, 
of e = 2.40 • 10~ 3 . Thus, both methods solve the prob- 
lem with satisfying accuracy. The regularized model is 
however found to be substantially more stable. To make 
this statement more quantitative, a series of simulations 
is run, on which the velocity (and thus the Mach num- 
ber) is kept constant at urj = 0.02v. For several chosen 
grid sizes N, the maximal Reynolds number Re max at 
which the simulation remains stable (i.e. delivers finite 
numerical values) is determined. Figure 3 shows that, 
although both methods exhibit a linear relationship be- 
tween Re max and N, the observed increase rate is 7.7 
times higher for the regularized method than for BGK. 

Finally, the capacity of the regularized model to rep- 
resent 3D flows has been explored in a preliminary study 
on direct numerical simulations (DNS) of a homoge- 
neous and isotropic turbulent flow. The system pos- 
sesses periodic boundaries and is driven by an external 
force that excites two wavenumbers in the limit of large 
wavelengths [15, 16]. It is known that the energy in- 
jected in such a system is mainly dissipated at the small- 
est scales, whose size is estimated by the so-called Kol- 
mogorov length If.. If these scales are not resolved with 
sufficient accuracy in the simulation, the system accumu- 
lates the energy and develops numerical instabilities. Our 
numerical simulations show that indeed, when the Kol- 
mogorov length is of the order of magnitude of a lattice 
site, Ik = 0.5 8r, with an average velocity u = 0.04 v, both 
BGK and the renormalized model exhibit a numerically 
stable flow. Furthermore, their statistical properties are 
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